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Abstract 

In this paper we describe a method to generate amorphous structures with arbitrary structural 
constraints. This method employs the Simulated Annealing algorithm to minimize a simple yet 
carefully tailored Cost Function (CF). The Cost Function is composed of two parts: a simple 
harmonic approximation for the energy-related terms and a cost that penalizes configurations 
that do not have atoms in the desired coordinations. Using this approach, we generated a set 
of amorphous carbon structures spawning nearly all the possible combinations of sp, sp^ and sp^ 
hybridizations. The bulk moduli of this set of amorphous carbons structures was calculated using 
Brenner's potential. The bulk modulus strongly depends on the mean coordination, following a 
power law behavior with an exponent v = 1.51 it 0.17. A modified Cost Function that segregates 
carbon with different hybridizations is also presented, and another set of structures was generated. 
With this new set of amorphous materials, the correlation between the bulk modulus and the 
mean coordination weakens. This method proposed can be easily modified to explore the effects 
on physical properties of the presence hydrogen, dangling bonds, and structural features such as 
carbon rings. 

PACS numbers: 61.43.Bn, 61.43.-j, 62.20.de 
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INTRODUCTION 



Carbon is an impressively versatile chemical element. As it is found in three distinct 
hybridizations, sp^, sp^ and sp, each with a well-defined local topology, this element can 
form a variety of different allotropes. Diamond's bulk hardness and graphite's laminar 
softness, for instance, can each be tracked down to carbon sp^ tetrahedron-like bonding and 
sp^ planar configuration. Although these two materials exhibit quite distinctive physical 
and chemical properties, they represent only a small fraction of possible carbon solids. 

Due to their many industrial applications, non-crystalline carbon materials have lately 
received attention. This class of material can be cheaply produced by Chemical Vapor 
Deposition (CVD) [1] and deposited over surfaces as thin, hard films, which exhibit good 
biocompatibility and chemical inertness j2|. In these films, the amount of sp^, sp^ and sp 
hybridized carbon, along with the presence or absence of hydrogen, directly infiuences the 
coating's stiffness [sl]. Considering the immense variety of amorphous carbon films that 
can be generated experimentally, many efforts were made to theoretically address how their 
properties vary according to changes in composition and local structure. 

One of the first models to describe non-crystalline materials was elaborated by Zachari- 
asen, who introduced the concept of continuous random network (CRN) to explain the 
atomic arrangement in Si02 glasses j4|. Despite primarily addressing Si02 structures, 
Zachariasen proposed that most oxide glasses could be considered as a random disposal 
of atoms in which there are neither bond defects nor long range crystallinity. Further ad- 
vances in this field were made possible by employing computers to generate continuous 
random networks. Among the most successful approaches is the one by Wooten, Winer and 
Weaire (WWW) [5]. Their ingenious bond-switching method consists of randomly swapping 
bonds from an originally 100% sp^ crystalline structure. In the WWW approach, one starts 
with a periodic diamond supercell, and follows a cycle of bond interchanges and geome- 
try relaxation until a fully tetrahedral amorphous carbon (ta-C) structure, also referred to 
as amorphous diamond (a-D), is obtained. This method is not only computationally fast 
and straightforward to implement, but also very successful in reproducing the experimental 
radial distribution function of amorphous sohds {g]. 

As computer performance improved, it became possible to generate a-C by employing 
Molecular Dynamics (MD). Using this strategy, it is possible to simulate either the quench- 
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ing of a melted carbon liquid or the process of film deposition, and amorphous carbon 
structures containing sp, sp^ and sp'^ hybridizations may be obtained jz-JJ]. Another ap- 
proach, analogous to MD, for the computer generation of amorphous carbon structures is 



the activation-relaxation technique, proposed by Barkema and Mousseau [13|. This method 
also allows one to obtain a CRN containing carbon atoms in different coordination, and it 
is an efficient means to overcome energy barriers between metastable minima. 

Many previous works have pointed out a general trend of a-C structures to become 



denser and stiffer with the increase of the mean atomic coordination 



These 



indings 



15| and by the percolation model of Phillips 



16| and 



are supported by experimental data 

Thorpe l3, 18]. The latter work may be seen as the theoretical basis that explains the 
strong dependence of elastic properties both on the mean coordination and on the presence 
of small rings, and it can also successfully explain CRNs' bulk modulus vanishing as the 
mean atomic coordination 1 approaches 2.4. 

Computer methods to generate a-C have recently included the use of increasingly sophis- 



ticated Hami 
atoms 
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tonians (such as ab initio) to more accurately simulate the dynamics of carbon 
22| |. Likewise, many of these algorithms focus on simulating the experimental 
process that originates a particular amorphous system. This may be seen as a top-down 
strategy, as for each new a-C material produced, it is necessary to perform an MD that 
resembles the corresponding experimental conditions. Only after that can the material's 
properties be estimated. For instance, if one were to answer the effect of a particular lo- 
cal structure on the amorphous carbon properties, various MDs might possibly have to be 
performed until a certain simulation yielded the desired microscopic structure. 

Later, with the discovery of new forms of a-C, the influence of some local scale features 



on the properties of amorphous systems was pointed out 
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24| . The intrinsic complexity 
of dealing with these features, such as the presence of rings, the proportion of sp, sp^ 
and sp^ hybridized atoms, and the size of clusters, may explain the difficulty to obtain an 
universal relationship among density, bulk modulus and mean coordination [231]. In order to 
efficiently study these aspects, it is important to develop a method to generate amorphous 
carbon structures, as neither WWW's algorithm nor MD can easily generate a-C with generic 
structural constraints. The former because it is limited to sp^ hybridized carbon, and the 
latter because it follows a top-down approach and offers no direct way to control the presence 
of those features. 
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In order to tackle these difficulties, we introduce a scheme for the computer generation of 
a-C structures which follows a bottom-up strategy. Relying on the algorithm of Simulated 
Annealing (SA) [25], we propose a Cost Function (CF) that aims not at simulating realistic 
atomic dynamics - as done in [26] - but rather at exploring several metastable configura- 
tions that meet some desired criteria. Our aim is to introduce something of a "theoretical 
workbench" to flexibly simulate a-C. One of our main interests is to employ this technique 
to calculate the bulk modulus' dependence on the fraction of sp, sp^ and sp^ hybridized 
carbon, and to perhaps find an amorphous structure more incompressible than diamond. 

The paper is organized as follows: First, in the computational strategy section, we detail 
the proposed Cost Function. Next, we present the results of simulations with 512 carbon 
atoms, using periodic boundary conditions, in which we calculate the bulk modulus of several 
amorphous structures spawning several possible combinations of sp^, sp^ and sp carbon 
and also compare the calculated radial distribution functions with the literature. We then 
show that we can modify our CF to force the creation of clusters of atoms with the same 
hybridization, and thus explore their effect on the bulk modulus. We conclude by discussing 
the advantages of the approach described in this paper, as well as possible extensions that 
can be proposed to deal with other problems. 

COMPUTATIONAL STRATEGY 

Following the idea of the SA technique, we have developed a fast and customizable CF 
which guides the exploration of different atomic configurations. Our Cost Function is com- 
posed of two parts, the ffist a computationally simple potential, and the second depending on 
the desired constraints for the CRN. The CF is required neither to be a continuous function 
nor to yield a realistic value for structures far from a metastable situation; its only require- 
ment is to have an arbitrary low value for geometrically stable configurations. Thus, the 
problem of finding a certain amorphous material that meets some arbitrary constraints can 
be transformed to the problem of finding a reasonably low minimum of this Cost Function 
using SA. 

Considering the potential part of the CF, a simple harmonic-like approximation was 
employed. As we wished our model to be as computationally fast as possible, carbon atoms 
were considered to be either bonded or non-bonded, with a cut-off distance Tc of 2.0 A. 
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This value is somewhat arbitrary, as long as it is smaller than the typical second-neighbor 
distance, but we found that shrinking it too much allows non-bonded atoms to remain too 
close. Thus, the potential term of the Cost Function is written as 



+v,j2[i~{u,-u,r] (1) 

Ui,Uj 

The first sum is over all bonds r^j and expresses the stretching energy relative to the 
equilibrium distance r*^i)c{j) between atoms i and j (having coordination c{i) and c(j) re- 
spectively). The second sum comprises all 6ijk angles having a common j center, where 6^^^^ 
denotes the equilibrium angle, which in our approximation depends only on the hybridiza- 
tion c(j) of atom j. The last sum considers only connected sp"^ centers, and it constitutes 
the torsional energy of having two sp^ planes, with normal vectors Uj and uj, nonparallel. 
In all summations, repeated terms are not counted. Also, the constants Vr, Va and Vt merely 
set the relative strength of the radial, angular and torsion terms for the total (pv 

The equilibrium quantities do not need to be found with great precision. Since our 
objective is simply to obtain a structure obeying a given set of constraints, there can be 
small geometrical distortions, all of which can be removed in a further relaxation by using a 
more realistic Hamiltonian. Thus, the following approximations were made: if two bonded 
atoms are both sp^, their equilibrium distance rl^ is the same found in diamond; for sp^-sp^ 
bonds, equilibrium interatomic distance rgg is that found in graphite and, in the case of sp-sp 



bonds, one takes the equilibrium distance of the triple bond in 2-butyne for [27|]. For a 
pair of bonded carbon atoms with different hybridizations, say, c' and c", the equilibrium 
distance r*,^„ is simply is the average (r*,^, +r*„^„)/2. Finally, if an atom has a coordination 
d > 4, the same value as for four-fold atoms is assumed (r*,^, = rl^). Likewise, the distance 
rli between singly-bonded atoms is the same as for sp-sp. 

The way we presented (pv alone will not yield any bondings, as linked atoms will always 
increase the system's energy. In order to correct this and to control the material hybridiza- 
tions, we introduce another term in the Cost Function. This term, here referred to as 
Coordination Cost, allows one to set how many atoms should be sp^, sp^ and sp: 
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0c = ^ ec'\nc' - n*/ 1 (2) 

c' 

The sum is over all possible c' coordinations, where ric' is the number of c'-coordinated 
atoms, and n*, is a parameter that sets how many atoms should be c'-coordinated. This 
way, each constant sets a cost for a configuration having a wrong number \nc' — of 
c'-coordinated centers. Clearly, we must have n*, = J2c' = N, where is the total 
number of atoms. Taking the absolute value of 12^' ~ n*, rather than squaring it has shown 
some advantages, such as the Cost Function exhibiting a sharper minimum. Moreover, it 
naturally makes the CF an extensive function, so that constants do not have to be altered 
as the number of atoms in the simulation box changes. 

Finally, the Cost Function $ is simply defined as a linear combination of the previous 
terms: 

$ = Xv(l)v + ^c4^c, (3) 

where Ay and Ac are constants. With the former definitions, bonded atoms are stable, 
provided that their linking decreases the number of atoms with wrong coordinations. By 
putting atoms in a cubic cell with periodic boundaries and setting how many should be 
s]?^, sp^ and sp hybridized (i.e., fixing 77,4, nl and a CRN can be obtained as the set 
of atomic positions which minimizes We expect our CF to possess many metastable 
minima, and we employed the Simulated Annealing (SA) algorithm to optimize 

Despite being SA a technique aimed at finding the global minimum of complex hypersur- 
faces, we are not necessarily interested in using this technique to its full extent. Consider, 
for instance, the generation of fully sp^ carbon CRNs. The global minimum of the CF in 
this case is crystalline diamond, which is clearly not the solution we are looking for. Ac- 
cordingly, we propose the use of the SA algorithm to find deep minima of the Cost Function 
that comply with a certain set of constraints. Whether or not this minimum is the global 
minimum of the Cost Function is not our concern. 

We devised the optimization algorithm the following way: Each step in the SA scheme 
constitutes randomly displacing one atom inside a cube of side length 0.4 A, changing the CF 
by A$. Even though system-wide movements could be implemented, individual movements 
have shown a great advantage, as $ only changes locally in this case. So, the recalculation 
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TABLE I: Parameters used in the three-regions anneahng scheme: initial (Pj) and final (Pf) 
acceptance rate and the fraction of the total steps assigned to each region. 



Region Steps P, Pf 

1 10% 95% 70% 

2 80% 70% 30% 

3 10% 30% 5% 



of (f)v, which is the most numerically intensive task in our CF, only has to be evaluated in 
a small radius around the displaced atom. This way, the computational complexity of A$ 
approximately does not scale with the system size. As in classical SA, each movement is 
accepted stochastically according to the weighing factor e"''^*, where (3 is the inverse of 
the fictitious temperature T. In addition to atomic displacements, the system periodically 
undergoes a random expansion or contraction, so that one does not need to set the final 
density a priori. Although A$ can not be evaluated locally in this case, the number of atomic 
displacements between full system scalings are such that these resizings do not compromise 
the algorithm speed. 

A key component to successfully finding a minimum for $ is the determination of the 
optimal annealing scheme 281]. Following Christoph and Hoffmann 291], we decreased T 
using a power law. In order to further control the annealing, we separated it in three 
regions, and following Johnson et al. [28], we fixed the initial and final acceptance rates 
instead of the temperatures, as the former proved to be less sensitive to changes in the CF. 
The initial and final acceptance rates and number of steps assigned to each annealing region 
were found by minimizing the width of the angular distribution function for a-D, and results 
are summarized in Table [B 

The constants in Eqs. ([I])-© were determined as follows. Medium-sized (256 atoms) 
mixed-coordination amorphous carbon networks were generated, and the resulting structures 
were further relaxed using Molecular Dynamics and Brenner's interatomic potential [sol 



using GULP 



311 ]. Parameters were thus varied to minimize the final structures' energies 



and coordination errors {i.e., the numbers of c'-coordinated atoms should be as close as 
possible to the defined number Uc'*)- The optimal constants are shown in Table HTl Also, 
we found Ay/Ac ~ 0.75 to be appropriate for most simulations. One should note that 
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TABLE II: Parameters for Eqs. ([H)-®. 



Constant 


Ay 


Ac 


Vr 




Vt 




Value 


1.0 


2.5 


5.0 


3.0 


1.5 




Constant 


eo 


ei 


£2 


£3 


£4 


" J vj — 


Value 


10.0 


5.0 


2.0 


1.5 


1.0 


10.0 lO-'' 


Constant 


Ml 


r* 
'22 


'33 


'44 


'^ii(i>4) 




Value (1) 


1.2 


1.2 


1.42 


1.54 


1.54 




Constant 


ei 




01 


ft* 

j(i>4) 






Value 


180° 


120° 


109.4° 


109.4° 







these constants could not be determined with great precision due to the large statistical 
fluctuations in the methodology. Fortunately, precise values are not sought, since small 
stresses in the CRNs could be eliminated a posteriori by a supplemental Molecular Dynamics 
using more sophisticated interatomic potentials, such as Brenner's potential 30| , or even ah 
initio calculations. 

As a next step, in order to validate our CF, we first generated a 64-atom 100% sp^ CRN. 
We used this CRN and a diamond structure as references to build 118 other CRNs, each of 
them constructed as a linear combination of the two reference structures. More specifically, 
denoting and the positions of the ith. atom of the diamond and the amorphous structure 
respectively, each interpolated CRN was defined according to Vi{u) = (1 — m) rf + ur^, where 
u is an interpolation parameter. 

For each CRN, we calculated the energy using our Cost Function and Brenner's potential, 
and plotted the results in Fig. [H It is clear that both models should yield high values 
for unstable materials {i.e., for u far from or 1), and that our potential should only 
partially reproduce the true energy surface. Nevertheless, our simplified Cost Function very 
much resembles the computationally more expensive Brenner's potential and, particularly, it 
agrees with it on the position of the two minima, which is the key feature for the Simulated 
Annealing scheme to identify a metastable CRN. 
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FIG. 1: (Color online). 



Comparison of the energy calculated using our Cost Function (top red 



curve) and Brenner's [30| potential (bottom blue) for a set of structures. Top curve has been 
displaced vertically for clarity. The left minimum at u = corresponds to a diamond structure, 
and the right one at n = 1 to a sp^ amorphous material. 

RESULTS AND DISCUSSION 



The procedure we are proposing for the generation of customized continuous random net- 
works was first apphed to map the bulk modulus' dependence on the fraction of sp, sp^ and 
sp^ hybridized carbon. To do so, we generated 45 amorphous structures, each containing 
512 atoms and having a different proportion of the possible atomic coordinations. These 
proportions were chosen so as to homogeneously cover a ternary graph mapping all possible 
coordinations. It took 2 x 10^ iterations (4 hours[425) for each SA simulation. After these 
structures were generated, they were submitted to Molecular Dynamics simulations using 
Brenner's potential to minimize non-physical features, such as distorted angles and dis- 
tances. The MD also ensures that each structure stays in a relatively low energy metastable 
configuration. 

Each MD simulation was carried out at a low temperature of 50 K in order to preserve 
the main features of the CRN generated by SA. The structures that were modified the most 
by the Molecular Dynamics process had roughly 50% sp/sp^ carbon atoms, but little or no 
sp"^ centers. In the worst case, a structure with 'z = 2.87 had 14.65% of its sp^ and 2.15% 
of its sp atoms turned into the sp'^ form. On the other hand, the most sp'^-rich structure, 
having final mean coordination z = 3.98, had less than 1.5% change in its hybridization due 
to the MD relaxation. 
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Due to the high computational cost required to perform a full MD, only the equilibration 
phase was performed. Each MD isobaric simulation was performed for 5 ps using a 0.1 fs time 
step. Afterwards, each structure was submitted to a full Hessian-driven geometry relaxation, 
so that the elastic moduli could be calculated. We could have also estimated the the elastic 
moduli directly using the volume fluctuations of the MD, but it would have required a much 
longer time period. Both the MD simulations and the Hessian-driven geometry relaxations 
were performed using GULP [sil]. Fig. |2] shows some CRNs generated, including a sp 
rich carbon network, which may be quite difficult to obtain using Molecular Dynamics 
under conventional approaches. One should note, however, that Brenner's potential does 
not includes van der Waals forces, which should be quite important in determining the 
geometries and elastic properties of these low-density sp carbon structures. 

As an example of the flexibility of our approach to generate a-C, we also proposed another 
term to the Cost Function to penalize binded atoms with different coordinations. This term 
was proposed to segregate carbon atoms having different hybridizations for main reasons. 
The first is purely theoretical, as, at least in principle, sp'^ centers embedded in a rigid sp'^ 
matrix should not infiuence the bulk modulus much, whereas if these threefold atoms form 
a segregate phase they might make the material considerably more compressible. This sp'^ 
segregation is also motivated by the existence of experimental carbon materials containing 



sp^-rich phases and graphite-like sp^-rich regions 32|]. The second reason is that it has 
been pointed out that the microstructure of hydrogenated a-C, mainly the size and shape 
of sp^ clusters, plays an important role in determining the electronic properties of these 
materials [33|. So, we proposed an heterogeneity term (pn, 



(pH = - Scii),c{j)) (4) 

where we are again summing over all possible bonds, and c{i) is the coordination of the 
center i, and 6 is the Kronecker delta function. 

With the addition of this term, the cost function becomes: 

$ = Xv(j)v + >^c(f)c + >^H(j)H (5) 

In our simulations, we used Xh = 1-5, which was the smaller number for which atoms 
with different coordinations would be visibly segregated. This by no means precluded that 
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(a) 




(c) 

FIG. 2: (Color online). Example of some generated structures. The following color scheme was 
used: sp^ atoms are shown in green, sp'^ ones in blue and sp in red. (a) A sp rich amorphous 
carbon network, (b) A mixed sp'^/sp^ structurl Vith Xh = 1.5. Note that the CRN segregates sp"^ 
and sp^ atoms in two phases due to the heterogeneity cost, (c) Another mixed sp^/sp^ structure 



a small number of atoms with other coordinations were found in the segregated phases. One 
may control how pure in terms of hybridization these regions are by varying Xh- Using 
this new $, we generated another set of 45 structures, later employing Brenner's potential 
as well. In order to distinguish from the other set of structures, we shall call the latter 
a-Cs (generated with Xh > 0) heterogeneous structures, and the former ones (with Xh = 0) 
homogeneous structures. Fig. [2] shows the effect of including (pn into the Cost Function by 
showing two CRNs generated with different values of but having similar amounts of sp^ 
and sp^ centers. 

The possibility of creating homogeneous and heterogeneous structures exemplifies the 
flexibility of our method to generate amorphous structures. By adding a simple and intuitive 
term to $, it is possible to generate CRNs with quite different characteristics. This same 
approach can be employed to generate CRNs having other microscopic features. For instance, 
a term may be added to increase the energy of a CRN if n-fold carbon rings are present. 
Also, the constant may have a non-zero value, and centers with only one bond may be 
readily mapped into hydrogen atoms or dangling bonds. 

One important question that can be readily answered using our algorithm is how the 
bulk modulus varies with the atomic coordinations. Although it has been pointed out that 
the bulk modulus should depend mainly on the mean coordination, no previous method 
exists to generate a-C with predetermined fractions of sp^, sp^ and sp carbon, so that this 
aspect has not yet been fully studied. We show in Fig. [3] the resulting bulk modulus as a 
function of the fraction of sp^, sp^ and sp carbon, for both the cases of homogeneous and 
heterogeneous structures [431]. We also plot the bulk modulus versus the mean coordination 
for the homogeneous case in Fig. |H The largest bulk modulus we found was 362 GPa, which 
is lower than that calculated for crystalline diamond using Brenner's potential (442 GPa, 
the same as the diamond's experimental bulk modulus 3J]). 

For the homogeneous case, the bulk modulus depended mainly on the mean atomic 
coordination, with a Spearman's rank correlation coefficient 351] p = 0.98 - supporting 



previous studies also pointing out this trend [12 
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isl. 



For heterogeneous networks, the 



dependence on the mean correlation diminished a little, with p = 0.9. However, considering 
only the region with mean coordination 2.5 < z < 3.5, both correlations drop to p = 
0.94 and p = 0.83, respectively. We explain the larger decrease of p for heterogeneous 
structures structures this way: since sp hybridized atoms form floppy 17| phases with null 
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B (GPa) 




FIG. 3: (Color online). Bulk modulus dependency on carbon hybridization. In each triangle, the 
lower left vertex represents a 100% sp structure (having a mean coordination z = 2), the top vertex 
a 100% sp'^ structure (with z = 3), and the lower right vertex a 100% sp^ material (with z = 4). 
Points lying on the same vertical line have the same mean coordination. Top: No constraint 
was imposed on the heterogeneity {Xh = 0). The bulk modulus varies little along any vertical 
line, suggesting that it may be well described by the mean coordination. Bottom: Heterogeneous 
structures generated with Xh = 1.5. The mean coordination does not dictate the bulk modulus as 
well as in the previous case, since B varies along vertical lines. 
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FIG. 4: (Color online). Variation of the bulk modulus as a function of the mean coordination. Light 
diamonds (dark circles) represent data from CRNs generated b y S A without (with) heterogeneity 



12l |. Inset: solid line representing 



cost. For comparison, crosses show tight-binding results from 
the power law fit to data for homogeneous CRNs. 

bulk modulus, some heterogeneous CRNs, such as one made of 50% sp^ and 50% sp carbon, 
will have a very small bulk modulus due to the large floppy region. Such small bulk modulus 
will not be observed in a 100% sp^ network, even though both structures have the same mean 
coordination. If we put \h = and let homogeneous structures form, the sp carbon will 
not segregate, but it will be incorporated between sp^ centers. Thus, there will be no large 
floppy regions. 

The bulk modulus is also plotted as a function of the mean coordination (Fig. HJ. Fol- 



lowing previous studies 
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18| . we fitted a power law to the bulk modulus data for the 



set of homogeneous CRNs, 



B{-z) = B,{z--Zfr (6) 

We found the phase transition from rigid to floppy states [l7| to occur at mean coor- 
dination Zf = 2.10 ± 0.11, with Bo = 140 ± 26 GPa and_u ^1^51 ± 0.17. These results. 



particularly the exponent, are close to those reported by [12|, [Ij, [18], as compared in Table 
mil The slight deviation for If can be explained by the size of the simulation cell: Even for 
relatively large cells containing 512 atoms, there is a chance that a non-floppy carbon chain 
of sp"^ or sp^ atoms will percolate the periodic cell. This was not observed by Mathioudkis 
et al. - whose results were extrapolated for mean coordination bellow 1 = 2.68 - nor by 



He and Thorpe [18|] and Djordjevic and Thorpe 



14| . because of a limitation of the bond 
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TABLE III: Comparison of the fitted parameters for Eq. 



Reference 



He and Thorpe [18] 2.4 1.5 ± 0.2 

Djordjevic and Thorpe [141 2.4 1.4 

Mathioudakis a/. [12j 2.33 1.5 ±0.1 

This Work 2.10 ±0.11 1.51 ±0.17 



depleting method which causes the simulation cell to collapse for small 

It must be pointed out that long-range effects may be taken into account after the amor- 
phous structures are generated. In our case, we employed the Brenner potential, which does 
not include such interactions, for the relaxation and calculation of the bulk moduli. It is 
reasonable to assume that the long sp chains are weakly binded by dispersive forces, so that 
the bulk modulus of floppy networks do not vanish completely. So, it is quite possible that 
using a potential model for the calculation of the elastic properties that includes van der 
Waals interactions would yield higher bulk moduli for low 1. 

After generating the CRNs, the set of 90 structures may be seen as the initial set of an 
expanding database that may be used, for instance, to extract structural information from 
experimental results. As an example of this application, we compared the calculated radial 
distribution functions (RDF) for the set of CRNs with results from the literature in Fig. |5l 
Using all our available structures, we searched for CRNs that would best reproduce the 
RDF of some experimental materials: one sputtered a-C 361] and one ta-C 371]. The first 
experimental a-C was prepared by rf sputtering, while the ta-C was grown using filtered 
cathodic arc. Using least-squares fitting, we found that a 88% sp^ and 12% sp homogeneous 
structure best reproduced the RDF of the sputtered a-C, while a heterogeneous 50% sp^ / sp^ 
structure best described the ta-C one. Furthermore, by adding an additional degree of 
freedom to scale the r variable, the structure that best described the ta-C one was CRN 
containing 80% sp^, 10% sp^ and 10% sp carbon atoms. 

The fit is not optimal: The mean coordinations of the experimental structures were 
reported to be 3.34 and 3.9, while the fitted mean coordination were 2.9 and 3.5 (3.74 if we 
add the additional degree of freedom). This error could be related to the potential used in 
the relaxation process, or to the number of iterations used to generate the structures, which 
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Radial Distance r (A) 



FIG. 5: (Color online). Reduced Radial Distribution Function G{r). Top curves: sputtered a- 
36 1 (dotted line) and best fit (100% sp^ structure, green solid line). Middle curves: ta-C [s^] 



C 



(dotted line) and best fit (heterogeneous 50% sp^ and 50% sp^ structure, red solid line). Bottom 
curves: the same experimental curve was used [s^], but the theoretical curve could also scale in 



the r axis. The yellow curve represents a CRN containing 
atoms. 



sp^ , 10% sp^ and 10% sp carbon 



in turn control their angular and bond distribution widths. However, even though the CRNs 
were not generated for this purpose, the comparison between experimental and theoretical 
RDFs was performed to show that the generated CRNs indeed present similarities with the 
experimental structures, to such an extent that they are able to reproduce experimental 
RDFs. It should be noted that it has never been in the scope of this work to present a 
method to extract structural information or create structures based on given experimental 



RDFs. There are specialized methods for this tasks, such as Reverse Monte Carlo [38| and 
Hybrid Reverse Monte Carlo 39|, which are much more efficient to extract information from 
experimental RDFs, but not to generate CRNs having particular coordination or structural 
constraints. 

Finally, we generated a set of structures to test the performance of the algorithm to 
generate 100% sp^ structures. CRNs having 64, 128, 256 and 512 atoms were generated, 
and two structures were created for a given number of atoms. After relaxation using the 
same strategy as before (molecular dynamics followed by Hessian-driven relaxation), the 
calculated angular widths ranged from 11.8° to 13.9°, with a mean value of 13.2°. We did 
not observe a significant correlation between the angular width and the number of atoms, 
and the angular distribution was relatively symmetric with an average of 109.07°. For 
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comparison, high-quality tetrahedral networks having an angular width of 9.19° have been 



assembled using a modified version of the WWW algorithm by Barkema and Mousseau |40|. 
Also, a reduction of the angular distribution width may be possible by increasing the number 
of steps during the SA or by extending the relaxation process after the structure is generated. 
So, although our strategy is not optimized for a-D as other methods, it is flexible enough to 
generate both a-D and a-C with various hybridizations and structural constraints. 

CONCLUSION 

In this paper, we described the computational creation of carbon CRNs employing the 
Simulated Annealing algorithm. We proposed a numerically simple Cost Function able to 
yield extremely different amorphous materials. As an example of the capabilities of our 
algorithm, we generated amorphous structures spawning nearly all possible combinations of 
sp'^, sp^ and sp hybridized centers, and then calculated their bulk moduli using Brenner's 



potential [30||. We were also able to easily modify our CF to create heterogeneous materials, 
in which atoms with the same hybridization tend to segregate. With the set of homogeneous 
structures, we observed a phase transition from floppy to rigid networks, and a power-law 
fltting of the bulk modulus dependency on mean atomic coordination was in close agreement 
with the literature. However, we noticed that the mean coordination z did not correlate with 
the bulk modulus of heterogeneous networks as well as it did for homogeneous ones. This 
indicates that the heterogeneity may play a very important role in dictating the elastic 
properties of a-C. 

The strategy we described is completely universal and customizable, and modiflcations 
can be easily made to include other chemical elements, such as hydrogen, and to control the 
presence of other features, such as rings and dangling bonds. Once CRNs with particular 
features are generated, their physical properties can be estimated using more sophisticated 
Hamiltonians, including ab initio calculations, whenever it is computationally feasible. Fur- 
ther extension of this approach to include microstructural constraints in t he p rocess of 



generating CRNs (such as those possibly responsible for ultrahigh hardness in 4l|]) depends 
only on the availability of suitable computational resources. 
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